function [f, g, theta, xi] = rcnl_gmmobj_demand( ...
    dparams,           ...
    s_jt,              ...
    xlinear,           ...
    xnonlin,           ...
    cdindex,           ...
    dfull,             ...
    IVD,               ...
    W,                 ...
    IVD_alt,           ...
    W_alt,             ...
    fe,                ...
    expmval_mat,       ...
    replace,           ...
    rho_transform,     ...
    tolerance          ...
)

Pi  = dparams(1:end-1);
rho = dparams(end);
if rho_transform
    rho = rho_transform * (1/(1+exp(-rho)));
end

[delta, ~] = rcnl_meanval(Pi, rho, s_jt, xnonlin, cdindex, dfull, expmval_mat, replace, tolerance);

xi = NaN;
theta = NaN(size(xlinear, 2), 1);
if max(isnan(delta)) == 1
    f = 1e+10;
    g = NaN(size(dparams));
else
    ddelta = rcnl_ddelta( ...
        Pi,           ...
        rho,          ...
        delta,        ...
        xnonlin,      ...
        cdindex,      ...
        dfull,        ...
        rho_transform ...
    );

    [theta, xi] = f_gmm(delta, xlinear, IVD_alt, W_alt, fe);
    J = rcnl_jacobian(ddelta, xlinear, IVD, W_alt, fe, IVD_alt);

    dMoments = xi' * IVD;
    f = (dMoments * W * dMoments')/2;
    g = (J * W * dMoments');
end
